LPS MIA accute all samples

knitr::opts_chunk$set( 
  warning=FALSE,
  message=FALSE,
  results = 'asis',
  error = FALSE,
  tidy = TRUE,
  fig.show = "hold")

LPS-MIA base analysis

library(DESeq2)
library(xlsx)
library(Hmisc)
library(dplyr)
library(viridis)
library(pheatmap)
library(variancePartition)
library(RColorBrewer)
library(pcaExplorer)
library(vsn)
library(EnhancedVolcano)
library(rapport)
library(factoextra)
library(broom)
library(readr)
options(check.names = FALSE)

Load data

dir <- file.path("D:", "R", "R workspace")
setwd(dir)

load("04222021_LPS orgs_base.RData")

Update functions

PCA.covar <- function(df, meta, heatmap.r = F, data.r = F, heatmap.p = T, type = "spearman") {
    l <- length(colnames(meta))
    l2 <- l + 1
    pca.cov <- prcomp(t(df))
    PC.covs <- pca.cov$x[, 1:l]
    PC_cov_correlation <- rcorr(as.matrix(PC.covs), as.matrix(meta), type = type)
    PC_variance.explained <- PC_cov_correlation$r[1:l, l2:length(rownames(PC_cov_correlation$r))]
    PC_cov.cor.p <- PC_cov_correlation$P[1:l, l2:length(rownames(PC_cov_correlation$r))]
    if (heatmap.r) 
        pheatmap(PC_variance.explained, cluster_rows = F, cluster_cols = F, show_colnames = T, 
            color = heatmap.color.code)
    if (heatmap.p) 
        pheatmap(PC_cov.cor.p, cluster_rows = F, cluster_cols = F, show_colnames = T, 
            color = colorRampPalette(c("red", "white"))(50), display_numbers = T)
    if (data.r) 
        return(PC.covs)
}

createDT <- function(DF, caption = "", scrollY = 500) {
    data <- DT::datatable(DF, caption = caption, extensions = "Buttons", options = list(dom = "Bfrtip", 
        buttons = c("copy", "csv", "excel", "pdf", "print"), scrollY = scrollY, scrollX = T, 
        scrollCollapse = T, paging = F, columnDefs = list(list(className = "dt-center", 
            targets = "_all"))))
    return(data)
}

First step QC

0. Subsetting and QC

# Subsetting a df to samples of interest (accute)
subset.df <- gfilt.df[colnames(gfilt.df) %in% rownames(filt.meta[filt.meta$`Accute/chronic` %in% 
    c("Accute"), ])]
subset.df <- subset.df[rowSums(subset.df) >= length(colnames(subset.df)), ]  #Also removes ERCC SIs

subset.meta <- filt.meta_clean[rownames(filt.meta_clean) %in% colnames(subset.df), 
    ]

# Sort columns
subset.meta <- subset.meta[sort(rownames(subset.meta)), ]
subset.df <- subset.df[sort(colnames(subset.df))]

all(rownames(subset.meta) == colnames(subset.df))

[1] TRUE

# Quick density QC 550 x 450
plot.density(subset.meta, to.color = "Stimulation", to.plot = "`Library size`")

# Quick covariate analysis
cov_for.cor <- cov_for.cor[sort(rownames(cov_for.cor)), ]
cov_for.cor.sub <- cov_for.cor[rownames(cov_for.cor) %in% colnames(subset.df), ]
cov_for.cor.sub <- cov_for.cor.sub[!colnames(cov_for.cor.sub) %in% c("Accute/chronic")]
covar.cor.sub <- rcorr(as.matrix(cov_for.cor.sub), type = "pearson")
covar.cor.sub_1 <- covar.cor.sub$r

# 500x450
pheatmap(covar.cor.sub_1, color = heatmap.color.code)

# D100 PCA-covariate analysis 500x500
D100.PCAcovar <- PCA.covar(subset.df, cov_for.cor.sub, data.r = T, type = "pearson")
# PCA_cov_cor_R(cov_for.cor.sub, subset.df)

Normalization

1. Performing normalization

# For correcting the data
subset.meta$"Log(library size)" <- scale(subset.meta$`Library size`)
subset.meta$`Scaled percent. MT` <- scale(subset.meta$"Percent. MT")

# Create a DESeq2 matrix
dds.complex <- DESeqDataSetFromMatrix(countData = as.matrix(subset.df), colData = data.frame(subset.meta), 
    design = ~Experiment + Log.library.size. + Percent..ERCC + Scaled.percent..MT + 
        Stimulation)


# version 2 based on dds object Estimate library size correction scaling factors
dds.complex <- estimateSizeFactors(dds.complex)

vsd.complex <- vst(dds.complex)

Post-normalization QC

2. Performing post-normalization QC:

# Doing QC on the vst-transformed values
meanSdPlot(assay(vsd.complex))

sampleDists <- dist(t(assay(vsd.complex)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- colnames(vsd.complex)
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)

annot_df <- data.frame(colData(dds.complex)$Stimulation, colData(dds.complex)$Cell.line)
colnames(annot_df) <- c("Stimulation", "Cell line")
rownames(annot_df) <- colnames(dds.complex)

annot_colors <- vector("list", length = 1)

annot_colors$Stimulation[["CTR"]] <- "#00AFBB"
annot_colors$Stimulation[["LPS"]] <- "#E7B800"

annot_colors$`Cell line`[["OH1.5"]] <- "#2E7281"
annot_colors$`Cell line`[["OH2.6"]] <- "black"
annot_colors$`Cell line`[["PRO1"]] <- "#50B8CF"

# Heatmap intersample-distance
pheatmap(sampleDistMatrix, clustering_distance_rows = sampleDists, clustering_distance_cols = sampleDists, 
    col = colors, annotation_row = annot_df, annotation_colors = annot_colors)

# IQR of samples
norm.df <- assay(vsd.complex)
sOutliers.norm <- IQRoutlier.detector(norm.df, barplot = T)

# PCA 550x450
pca.norm <- plotPCA.custom(as.matrix(norm.df), intgroup = c("Cell line", "Stimulation", 
    "Experiment"), ntop = 50000, returnData = TRUE, metadata = subset.meta)
PoV.norm <- round(100 * attr(pca.norm, "percentVar"))
PCAplot(pca.norm, Condition = "Stimulation", PoV.df = PoV.norm)
PCAplot(pca.norm, "Cell line", Shape = "Experiment", PoV.df = PoV.norm)

Data correction

3. Correcting the data

# Remove batch effects Create a meta df that corresponds to the sorted columns of
# the vsd assay but has covariates as numeric factors
covs_forremBatch <- subset.meta
covs_forremBatch$Stimulation <- as.factor(as.numeric(covs_forremBatch$Stimulation))
covs_forremBatch$Experiment <- as.factor(as.numeric(covs_forremBatch$Experiment))
covs_forremBatch$`Cell line` <- as.factor(as.numeric(covs_forremBatch$`Cell line`))
covs_forremBatch$Batch <- as.factor(as.numeric(covs_forremBatch$Batch))
covs_forremBatch$Replicate <- as.factor(as.numeric(as.factor(covs_forremBatch$Replicate)))
covs_forremBatch <- covs_forremBatch[!colnames(covs_forremBatch) %in% c("Duplicate", 
    "Accute/chronic")]

batch.rem <- removeBatchEffect(as.matrix(assay(vsd.complex)), batch = covs_forremBatch$Experiment, 
    covariates = as.matrix(cbind(covs_forremBatch$`Log(library size)`, covs_forremBatch$`Percent. ERCC`, 
        covs_forremBatch$`Scaled percent. MT`)), design = model.matrix(~covs_forremBatch$Stimulation))
batch.rem2 <- removeBatchEffect(as.matrix(assay(vsd.complex)), batch = covs_forremBatch$Experiment, 
    batch2 = covs_forremBatch$Replicate, covariates = as.matrix(cbind(covs_forremBatch$`Log(library size)`, 
        covs_forremBatch$`Percent. ERCC`, covs_forremBatch$`Scaled percent. MT`)), 
    design = model.matrix(~covs_forremBatch$Stimulation))

Coefficients not estimable: batch24 batch25 batch26 batch27

Post-correction QC

4. Post-correction QC

# ====PC loadings for covariates====# 500x500
covs_forremBatch <- covs_forremBatch[colnames(covs_forremBatch) %in% colnames(cov_for.cor.sub)]
D100.PCAcovar.postcorrection <- PCA.covar(batch.rem, covs_forremBatch, data.r = T, 
    type = "pearson", heatmap.r = T)
# PCA_cov_cor_R(covs_forremBatch, batch.rem)

# PCA (550 x 450)
pca.cor <- plotPCA.custom(as.matrix(batch.rem), intgroup = c("Stimulation", "Cell line", 
    "Experiment"), ntop = 50000, returnData = TRUE, metadata = subset.meta, pc.1 = 1, 
    pc.2 = 2)
PoV.cor <- round(100 * attr(pca.cor, "percentVar"))

# 600x450
PCAplot(pca.cor, "Stimulation", PoV.df = PoV.cor)
PCAplot(pca.cor, "Cell line", Shape = "Experiment", PoV.df = PoV.cor)

Something that might make you sad: this is how the data looks when we can remove duplicate correlation:

pca.cor2 <- plotPCA.custom(as.matrix(batch.rem2), intgroup = c("Stimulation", "Cell line", 
    "Experiment"), ntop = 50000, returnData = TRUE, metadata = subset.meta, pc.1 = 1, 
    pc.2 = 2)
PoV.cor2 <- round(100 * attr(pca.cor2, "percentVar"))
PCAplot(pca.cor2, "Stimulation", PoV.df = PoV.cor2)
PCAplot(pca.cor2, "Cell line", Shape = "Experiment", PoV.df = PoV.cor2)

# Inter sample distances post correction
sampleDists.postcor <- dist(t(batch.rem))
sampleDistMatrix.postcor <- as.matrix(sampleDists.postcor)
rownames(sampleDistMatrix.postcor) <- colnames(vsd.complex)
colnames(sampleDistMatrix.postcor) <- NULL


# Heatmap intersample-distance 650x500
pheatmap(sampleDistMatrix.postcor, show_rownames = F, clustering_distance_rows = sampleDists.postcor, 
    clustering_distance_cols = sampleDists.postcor, col = colors, annotation_row = annot_df, 
    annotation_colors = annot_colors)

DESeq2 anaylsis

5. Performing DESeq2

dds.complex <- DESeq(dds.complex)

res.complex <- results(dds.complex, name = "Stimulation_LPS_vs_CTR")

summary(res.complex)

out of 17645 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 0, 0% LFC < 0 (down) : 0, 0% outliers [1] : 0, 0% low counts [2] : 0, 0% (mean count < 1) [1] see ‘cooksCutoff’ argument of ?results [2] see ‘independentFiltering’ argument of ?results

Interactive results table

(5.1) Results

# DT::renderDT(data.frame(res.complex), 'OH1.5A',scrollY=1000)

createDT(data.frame(res.complex), "Accute", scrollY = 1000)

Plotting results

6. Plotting results

#Volcano plot
#700x750
EnhancedVolcano(res.complex,
                lab = rownames(res.complex),
                x = 'log2FoldChange',
                y = 'padj',
                pCutoff = 0.05,
                FCcutoff = 2,
                labSize = 6,
                
                legendPosition = 'right')

# Grab significant downreg & highest lfc genes
if (length(res.complex[is.na(res.complex$padj), ]$padj) != 0) res.complex[is.na(res.complex$padj), 
    ]$padj <- 1

downreg.genes <- rownames(res.complex[res.complex$log2FoldChange < -2 & res.complex$padj, 
    ])
upreg.genes <- rownames(res.complex[res.complex$log2FoldChange > 2 & res.complex$padj, 
    ])


# Plot the genes heatmap based on sorted samples for Stimulation:
fact.to.sort <- as.numeric(annot_df$`Cell line`)
sorted <- batch.rem[, order(fact.to.sort, decreasing = F)]
sub <- sorted[rownames(sorted) %in% c(upreg.genes, downreg.genes), ]
# 550x1300
pheatmap(sub, cluster_rows = T, cluster_cols = F, annotation_legend = T, show_colnames = F, 
    annotation_col = annot_df, color = heatmap.color.code, scale = "row", annotation_colors = annot_colors)

LPS response & microglia genes

(6.1) LPS response genes

LPS.genes <- sort(c("IFT1", "IFT2", "IFT3", "CXCL8", "CX3CR1", "IL10", "IL1B", "TNF", 
    "IL6", "TMEM119", "AIF1", "TLR4", "IL12A", "IL23A", "CD40", "P2RY12", "GFAP"))
LPS.sub <- sorted[rownames(sorted) %in% c(LPS.genes), ]

pheatmap(LPS.sub, cluster_rows = T, cluster_cols = F, annotation_legend = T, show_colnames = F, 
    annotation_col = annot_df, color = heatmap.color.code, scale = "row", annotation_colors = annot_colors)

Patir et al. microglia markers

(6.2) Patir et al. marker genes

dir <- file.path("C:", "Users", "rapha", "Dokumente", "UM FN", "Internship", "Internship_Start", 
    "Project 4 - Synthesis", "Analysis", "Genelists")

setwd(dir)

coresign <- read_lines("PatirEtAl_coreMicrogliaGenes.csv")
genesOI <- batch.rem[rownames(batch.rem) %in% coresign, ]

pheatmap(genesOI, cluster_rows = F, show_rownames = F, cluster_cols = T, annotation_legend = T, 
    show_colnames = F, annotation_col = annot_df, fontsize = 6, color = heatmap.color.code, 
    scale = "row", annotation_colors = annot_colors)

pheatmap(genesOI, cluster_rows = F, show_rownames = F, cluster_cols = F, annotation_legend = T, 
    show_colnames = F, annotation_col = annot_df, color = heatmap.color.code, scale = "row", 
    annotation_colors = annot_colors)

Patir et al. microglia markers

7. Composite microglia marker score:

Patir_score <- prcomp(t(genesOI))

# identical(colnames(genesOI), rownames(subset.meta))


# rcorr(as.matrix(Patir_score$x), as.matrix(as.numeric(subset.meta$Stimulation)))

PCA.covar(genesOI, covs_forremBatch, data.r = T, type = "pearson", heatmap.r = T)

                    PC1        PC2         PC3        PC4         PC5

OH1-5B4D38CTR1 -1.74636258 2.8908017 1.02848364 2.5834555 -2.19440657 OH1-5B4D38CTR2 -1.00208645 0.6358254 -0.01398362 -0.4047980 2.03101594 OH1-5B4D38CTR3 -3.41508531 2.2693310 1.99621374 -1.1706207 -0.43971717 OH1-5B4D38LPS1 3.19382859 -2.1848698 -0.22014013 -0.5797534 0.88232874 OH1-5B4D38LPS2 0.54966486 -1.9875982 -4.72708866 1.6562309 -0.38398223 OH1-5B4D38LPS3 2.42004089 -1.6234901 1.93651503 -2.0845143 0.10476129 OH1-5B5D38CTR_1 -2.62972927 0.5794867 1.14752562 0.5710132 0.09500758 OH1-5B5D38CTR_2 0.75762124 1.3138442 -2.50090371 0.9778932 0.41035087 OH1-5B5D38CTR_3 1.67752849 1.3791259 0.29815748 -1.4062183 2.09530128 OH1-5B5D38LPS_1 2.07546124 -0.4512134 -0.12650997 1.0443274 0.12716672 OH1-5B5D38LPS_2 -0.53319395 -0.9457336 0.03402299 -1.6296096 -0.70720060 OH1-5B5D38LPS_3 -1.34768775 -1.8755098 1.14770758 0.4425942 -2.02062584 OH2-6B5D38CTR_1 10.17233506 2.0764612 0.41221900 -0.1750392 -0.92630455 OH2-6B5D38CTR_2 -4.01474252 1.2175377 -1.18427120 -1.0939976 1.02699449 OH2-6B5D38CTR_3 -3.19550551 1.3053340 -2.41939479 -1.4331930 -0.27125948 OH2-6B5D38LPS_1 0.04033274 -0.3485140 0.63748910 0.5721122 -0.94188988 OH2-6B5D38LPS_2 -2.45517285 -2.0958791 0.73582794 0.2377475 0.47652628 OH2-6B5D38LPS_3 -0.54724693 -2.1549398 1.81812994 1.8923702 0.63593315 PRO1B2CTR1D37 0.32692134 0.3375157 -0.67431100 -0.6502875 -0.77925089 PRO1B2CTR2D37 0.40579229 1.0141065 0.62233632 0.9588242 1.06056554 PRO1B2CTR3D37 0.61552080 1.1320469 -0.98752565 -1.0191720 -0.69814332 PRO1B2LPSD37_1 0.10321043 0.1890856 1.05778546 2.6107954 2.04463051 PRO1B2LPSD37_2 -0.71134949 -1.2447647 0.13640154 -1.6037324 -0.90367939 PRO1B2LPSD37_3 -0.74009537 -1.4279899 -0.15468666 -0.2964278 -0.72412245 PC6 PC7 PC8 OH1-5B4D38CTR1 0.21480100 -1.91113616 1.1044137 OH1-5B4D38CTR2 1.98553343 1.06121083 1.0745311 OH1-5B4D38CTR3 -0.25665357 2.36892821 -1.4395639 OH1-5B4D38LPS1 -2.02010049 -1.49275215 0.2772924 OH1-5B4D38LPS2 -0.02764314 0.48593103 -0.8854534 OH1-5B4D38LPS3 0.10406277 -0.51218176 -0.1312199 OH1-5B5D38CTR_1 -0.42263670 -0.08905834 -0.5954873 OH1-5B5D38CTR_2 0.94704612 0.50107298 0.4288915 OH1-5B5D38CTR_3 1.34975022 -2.15171211 -0.9893739 OH1-5B5D38LPS_1 0.53255664 0.77311571 0.1877774 OH1-5B5D38LPS_2 -1.99598448 1.42034357 2.2008616 OH1-5B5D38LPS_3 -0.41073179 -0.45376181 -1.2326692 OH2-6B5D38CTR_1 -0.22919664 1.10087039 -0.3718179 OH2-6B5D38CTR_2 -1.46489700 -0.64857151 -0.8236197 OH2-6B5D38CTR_3 -0.85455874 -0.43438780 0.2496625 OH2-6B5D38LPS_1 0.50967516 -0.19661800 0.8486596 OH2-6B5D38LPS_2 1.59006809 0.75444472 1.1424221 OH2-6B5D38LPS_3 0.44890913 -0.57573780 -1.0453066 PRO1B2CTR1D37 1.42461907 -0.23824278 0.5731002 PRO1B2CTR2D37 -1.23915952 -0.12769537 0.6874805 PRO1B2CTR3D37 -0.01835041 -0.54801741 -0.9117108 PRO1B2LPSD37_1 -1.61494444 0.72936799 0.2050794 PRO1B2LPSD37_2 0.77958128 -0.83222264 1.2139144 PRO1B2LPSD37_3 0.66825402 1.01681021 -1.7678636

PCA_cov_cor_R(covs_forremBatch, genesOI)

(7.1) Look at the PCA:

pca.mic <- plotPCA.custom(as.matrix(genesOI), intgroup = c("Stimulation", "Cell line", 
    "Experiment"), ntop = 50000, returnData = TRUE, metadata = subset.meta, pc.1 = 1, 
    pc.2 = 2)
PoV.mic <- round(100 * attr(pca.mic, "percentVar"))

PCAplot(pca.mic, "Stimulation", Shape = "Experiment", PoV.df = PoV.mic)

LPS response genes

8. Composite LPS response score:

LPS_score <- prcomp(t(LPS.sub))

# identical(colnames(LPS.sub), rownames(subset.meta))

# rcorr(as.matrix(LPS_score$x), as.matrix(as.numeric(subset.meta$Stimulation)))

PCA.covar(LPS.sub, covs_forremBatch, data.r = T, type = "pearson", heatmap.r = T)

                    PC1           PC2           PC3          PC4

OH1-5B4D38CTR1 1.55468250 -8.509066e-01 0.2200962835 0.644432419 OH1-5B4D38CTR2 0.74812227 7.748052e-01 0.2731713915 -0.388676504 OH1-5B4D38CTR3 1.64861884 2.307332e-01 -0.4379126279 0.077374434 OH1-5B4D38LPS1 -2.05432042 -5.194201e-01 0.0493785410 0.530878039 OH1-5B4D38LPS2 -2.00689777 -7.866440e-02 -0.4864094472 -0.716735003 OH1-5B4D38LPS3 0.10979458 4.434526e-01 0.3816758592 -0.147273384 OH1-5B5D38CTR_1 0.36203233 1.658608e-01 -0.4407139455 -0.099854301 OH1-5B5D38CTR_2 0.34246086 -6.347020e-01 -0.1715099740 0.005502484 OH1-5B5D38CTR_3 0.68271969 -4.371753e-01 0.9409558116 -0.430837597 OH1-5B5D38LPS_1 -0.73884468 3.413154e-02 0.3551640008 -0.082803648 OH1-5B5D38LPS_2 -0.73721914 6.912883e-01 -0.2504483797 0.095995761 OH1-5B5D38LPS_3 0.08885094 1.805967e-01 -0.4334475131 0.511997301 OH2-6B5D38CTR_1 -0.19871539 -4.323615e-01 1.0919632696 -0.137609629 OH2-6B5D38CTR_2 0.23770593 -4.683400e-01 -1.0401435226 0.059818441 OH2-6B5D38CTR_3 0.80368881 7.173346e-02 -0.4240476855 -0.783784105 OH2-6B5D38LPS_1 -0.51985909 6.578935e-01 0.5110810152 0.286818295 OH2-6B5D38LPS_2 -0.03137982 9.208090e-01 -0.2884259429 0.327056587 OH2-6B5D38LPS_3 -0.29144044 -7.497344e-01 0.1495728663 0.247700411 PRO1B2CTR1D37 0.02209927 5.255697e-01 0.5177850650 -0.269422233 PRO1B2CTR2D37 -0.21320435 -4.745589e-01 0.0002228491 0.001116586 PRO1B2CTR3D37 -0.03461903 -2.029023e-05 0.1518442785 0.686809971 PRO1B2LPSD37_1 0.07053913 -1.197639e+00 -0.3109711644 -0.504500690 PRO1B2LPSD37_2 0.23895682 7.799898e-01 0.1299434992 0.046036331 PRO1B2LPSD37_3 -0.08377184 3.666590e-01 -0.4888245274 0.039960035 PC5 PC6 PC7 PC8 OH1-5B4D38CTR1 0.06530352 0.589759011 -0.096116485 0.318414999 OH1-5B4D38CTR2 -0.34020462 -0.239409304 0.204750455 -0.219434232 OH1-5B4D38CTR3 0.23193953 -0.369539165 0.077466053 -0.369406890 OH1-5B4D38LPS1 0.07956146 -0.026517595 0.131620459 -0.078818005 OH1-5B4D38LPS2 -0.07548374 0.019550748 -0.330250024 0.189284851 OH1-5B4D38LPS3 0.03888386 0.026156305 0.012529542 0.159959277 OH1-5B5D38CTR_1 -0.26271326 -0.004211853 -0.359999938 0.209426238 OH1-5B5D38CTR_2 0.04058866 0.333723075 -0.006759743 -0.152903934 OH1-5B5D38CTR_3 -0.37911028 -0.197952352 -0.044080806 0.117434410 OH1-5B5D38LPS_1 0.17363790 0.084809491 0.232148431 -0.200006629 OH1-5B5D38LPS_2 0.44379713 -0.024740063 0.203899828 -0.018017546 OH1-5B5D38LPS_3 -0.01620016 -0.191628298 -0.025207772 0.044067460 OH2-6B5D38CTR_1 0.47491100 0.047839174 -0.393969929 -0.363118366 OH2-6B5D38CTR_2 0.51752801 0.033429169 -0.111469115 -0.142540604 OH2-6B5D38CTR_3 0.07118162 -0.110000696 -0.134063150 0.275588745 OH2-6B5D38LPS_1 -0.54446154 0.332350731 0.048679161 0.001208301 OH2-6B5D38LPS_2 -0.39688074 0.217088703 0.111195712 -0.079313477 OH2-6B5D38LPS_3 -0.12227836 -0.520707080 0.479627320 0.308175401 PRO1B2CTR1D37 0.50045590 0.170479557 -0.097363788 -0.014719874 PRO1B2CTR2D37 0.01822542 -0.131622565 -0.138240993 0.199186210 PRO1B2CTR3D37 -0.21257723 -0.593986294 -0.476261588 -0.073491228 PRO1B2LPSD37_1 -0.37391206 0.147361642 0.483277987 -0.269211806 PRO1B2LPSD37_2 0.54464456 0.081617562 0.394699435 0.342763697 PRO1B2LPSD37_3 -0.47683659 0.326150098 -0.166111052 -0.184526999

# PCA_cov_cor_R(covs_forremBatch,LPS.sub)

pca.LPS <- plotPCA.custom(as.matrix(LPS.sub), intgroup = c("Stimulation", "Cell line", 
    "Experiment", "Replicate"), ntop = 50000, returnData = TRUE, metadata = subset.meta, 
    pc.1 = 1, pc.2 = 2)
PoV.LPS <- round(100 * attr(pca.LPS, "percentVar"))

PCAplot(pca.LPS, "Stimulation", Shape = "Experiment", PoV.df = PoV.LPS)

Association between LPS response and microglia signatures

9. How strongly are the microglia signature and the LPS response correlated?

LPS.mic <- data.frame(cbind(LPS = scale(colMeans(LPS.sub)), `Microglia signature` = scale(colMeans(genesOI))))
colnames(LPS.mic) <- c("LPS", "Microglia signature")
library(ggpmisc)
formula <- y ~ x
ggplot(LPS.mic, aes_string(y = as.matrix(LPS.mic["LPS"]), x = as.matrix(LPS.mic["Microglia signature"]))) + 
    geom_point() + scale_colour_gradient(low = "gray90", high = "black") + stat_smooth(method = lm) + 
    theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method = lm, 
    linetype = "dashed") + labs(y = "Cumulative LPS response score", x = "Microglia signature score") + 
    stat_poly_eq(aes(label = c(paste(..adj.rr.label..))), label.x = "right", label.y = 0.37, 
        formula = formula, parse = TRUE, size = 3) + # geom_cor(method='pearson')+
stat_fit_glance(method = "lm", method.args = list(formula = formula), geom = "text", 
    aes(label = paste("P-value = ", signif(..p.value.., digits = 3), sep = "")), 
    label.x = "right", label.y = 0.3, size = 3)